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ABSTRACT 

We investigated the underlying architecture of planetary systems by deriving the distribution of planet multi- 
plicity (number of planets) and the distribution of orbital inclinations based on the sample of planet candidates 
discovered by the Kepler mission. The scope of our study included solar-like stars and planets with orbital peri- 
ods less than 200 days and with radii between 1 .5 and 30 Earth radii, and was based on Kepler planet candidates 
detected during Quarters 1 through 6. We created models of planetary systems with different distributions of 
planet multiplicity and inclinations, simulated observations of these systems by Kepler, and compared the prop- 
erties of the transits of detectable objects to actual Kepler planet detections. Specifically, we compared with 
both the Kepler sample's transit numbers and normalized transit duration ratios in order to determine each 
model's goodness-of-fit. We did not include any constraints from radial velocity surveys. Based on our best-fit 
models, 75-80% of planetary systems have 1 or 2 planets with orbital periods less than 200 days. In addition, 
over 85% of planets have orbital inclinations less than 3 degrees (relative to a common reference plane). This 
high degree of coplanarity is comparable to that seen in our Solar System, with the exception of Mercury. These 
results have implications for planet formation and evolution theories. Low inclinations are consistent with plan- 
ets forming in a protoplanetary disk, followed by a relatively calm evolutionary history without significant and 
lasting perturbations from other bodies capable of increasing inclinations. 

Subject headings: methods: statistical - planetary systems - planets and satellites: general - planets and satel- 
lites: detection 



1, INTRODUCTION 

Knowledge of the architecture of planetary systems can 
provide important constraints and insights into theories of 
planet formation and evolution. For instance, planetary sys- 
tems with a high degree of coplanarity or alignment, such as 
the Solar System, are consistent with the standard formation 
model of planets forming in a protoplanetary disk. Addition- 
ally, planetary systems with inclined or misaligned orbits can 
be indicative of past events that increased eccentricities and 
inclinations (e.g., Kozai oscillations by outlying perturbers, 
resonant encounters between planets, or planet-planet scatter- 
ing). Consequently, information on planetary multiplicity and 
the distribution of inclinations can reveal fundamental events 
in the lives of planetary systems as well as test theories of 
planet formation and evolution. 

Because of observational biases, it may be difficult to re- 
liably assess the underlying multiplicity of planets and their 
inclination distributions. Planets may go undetected if they 
have masses/radii below detection limits or if they are non- 
transiting (for the case of transit surveys). As for inclina- 
tions, mutual (or relative) inclinations between planets in the 
same system have only been measured in a limited number 
of cases or indirectly constrained. Exam ples include the plan- 
ets or biting the pulsar PSR B 1257+12 dKonacki & Wolszczanl 
2003), the planets GJ 876 b and c using radial velocity 
and/o r astrometric d ata dBean & Se ifahrt 20091 ICorreia et al.l 
120 lOt fealuev 2011), v And planets c and d using astromet- 
ric and radial velocity data (iMcArfhur et al.l l2010h . and us- 
ing transit timing/duratio n analyses to const rain inclinations 
in the Kepler-9 system dHolman et al.l l2010)), Kepler-10 sys- 
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tem (IBatalha et al.ll201 U iFressin et al.ll20TTl) , and Kepler- 11 
svstern lLissauer et al.ll201 lalh 

However, there now exists substantial knowledge about a 
large sample of planetary systems that can be used to statis- 
tically determine the underlying number of planets in each 
system and their relative orbital inclinations. In particular, 
the haul of planetary candidates catal oged by the Kepler mi s- 
sion (current coun t is over 2,300; Boruc kietaLlbo 11 alibi 
IBatalha et al.l2012l) provides a large trove of systems to study. 
These Kepler detections are termed planetary candidates, be- 
cause they have not been formally confirmed as real planets 
by independent observations such as radial velocity detections 
or successful elimination of false positive scenarios. We will 
refer to Kepler planetary candidates interchangeably as plan- 
etary "candidates" or "planets." 

In the present study, we use the latest Kepler catalog of 
planetary candidates (id e ntified in Quarters 1 through 6) re- 
leased by Batal ha et al.1 (120121) to constrain the underlying 
multiplicity and inclinations of planetary systems discovered 
by Kepler. By "underlying," we are referring to our esti- 
mate of the unbiased, true distribution. In our analysis, we 
create thousands of model populations, each obeying dif- 
ferent underlying multiplicity and inclination distributions. 
These planetary systems are subject to artificial observations 
to determine which planets could be transiting and detectable 
by Kepler. The detectable planets and their properties are 
compared to the observed sample to determine each model's 
goodness-of-fit to the data. Our results show that the underly- 
ing architecture of planetary systems is typically thin (low rel- 
ative inclinations) with few planets with orbital periods under 
200 days. Similarly, other statistical studies have been per- 
formed to constrain multiplicity and/or i nclinations using dif- 
ferent methodologies or data sets (e.g.. iLissauer et al.ll201 lbt 
iTremaine & Dong||2012t iFigueira et alj|2012t iFabrvckv et al l 



120 1 2t Uohansen e t al. 2012). Our analysis improves on previ- 
ous work by including all available quarters, extending to 200- 
day periods, and fitting models to observables such as normal- 
ized transit duration ratios that contain information on mutual 
orbital inclinations; these improvements lend to a deeper in- 
vestigation of the intrinsic distributions of planetary systems. 

This paper is organized as follows. Section [2] provides de- 
tails on how we created model populations including choice 
of stellar and planetary properties. Section [3] describes how 
we determined detectable, transiting planets from our simu- 
lated populations as well as our methods for determining each 
model's goodness-of-fit compared to observations. Section @] 
presents our main results with analysis of our best-fit models. 
Section [5] includes comparisons of our results with the Solar 
System and with previous work. Lastly, Section[6]summarizes 
the main conclusions of our study. 

2. CREATING MODEL POPULATIONS 

This section describes our methods for creating model pop- 
ulations. Each model population consists of approximately 
10 6 simulated planetary systems from which we determine 
planets detectable by Kepler. In total, we created thousands 
of model populations obeying different underlying planet (1) 
multiplicity and (2) inclination distributions. 

2.1. Stellar Properties 

Each model population's planetary systems have host stars 
with parameters drawn from a distribution of stellar proper- 
ties. For each simulated system, we randomly drew a star 
from a subset of t he quarterly KICQ (Kepler Input Catalog; 
iBrown et alj|201 lb CDPP (Combined Differential Photomet- 
ric Precision) lists. There is a separate list per quarter, because 
each file only lists targets observed that quarter and their cor- 
responding CDPP noise levels over 3, 6, and 12 hours. These 
quarterly lists also include basic KIC parameters (unchanging 
in the quarterly lists) for each observed star such as its Ke- 
pler magnitude K p , effective temperature T e ff, surface gravity 
parameter log(g), and radius 

We created a subset of the KIC CDPP lists in the following 
manner. First, we collected the KIC CDPP lists for Quarters 
1 to 6 (Ql-6). Not all stars are observed each quarter and we 
counted a total of 189,998 unique stars observed at some point 
during Ql-6. Next, we filtered this list of unique stars such 
that we only included stars observed under the "exoplanet" 
program, have nonzero values for and CDPP entries, and 
obey the following restrictions: 

4100 K < r eff < 6100 K, 

4.0 < log(g [cm s -2 ]) < 4.9, (1) 
K p < 15 mag. 

These stellar cuts given above in Equation (Q~|) allowed 
us to consider the brighter main-sequence stars with well- 
cha racterized propertie s in the KIC (consistent with those 
in iHoward et alj|2012h . With such filtering, we were left 
with a KIC subset of 59,224 remaining unique stars, which 
is roughly one-third of the original list of unique stars. For 
each star in this KIC subset, we calculated median CDPP val- 
ues over 3, 6, and 12 hours from all the values available in the 

1 The KIC is available for download from the 
Mikulski Archi ve for Space Telescopes (MAST) at 

|http : / / archive . stsci . edu/kepler/kic . html| 



Ql-6 lists. We also used each star's surface gravity parame- 
ter log(g) and R* values, which are common across the Ql-6 
KIC CDPP lists, to calculate each star's mass. Finally, we 
recorded each star's observing history over Ql-6 as indicated 
by the quarterly KIC CDPP lists. 

Our subset of the KIC with 59,224 stars was the catalog 
from which we randomly picked stellar hosts (with their cor- 
responding properties) for each simulated planetary system. 
In the next subsections, we describe how we chose physical 
and orbital parameters for the planets in each system. 

2.2. Planet Period and Radius Distribution 

We assigned values for the period and radius of each simu- 
lated planet by drawing from debiased distributions. To create 
such debiased distributions, we started with the observed sam- 
ple of Kepler candidate planets (or Kepler Objects of Interest; 
KOIs) and converted it to a debiased sample of planets using 
detection efficiencies. First we will discuss how we filtered 
the KOI catalog, and then we will discuss the calculation and 
application of detection efficiencies. 

Our observed sample of KOIs is based on the February 
2012 release (Ql-6) by iBatalha et alJ (120121) consisting of 
2,321 planetary candidates. We only considered a subset of 
this KOI catalog by applying the following filters. First, we 
removed any planet candidates with negative values for the 
following entries: orbital period, impact parameter, transit 
duration, semi-major axis, and signal-to-noise ratio (SNR). 
Negative values are given f or the period when a candidate has 
only one observed transit (Batalha et al.ll2012h . and we sus- 
pect that in general, negative values are assigned to quantities 
that have not been reliably measured. This cut left us with 
2,298 remaining candidates. Next, we filtered out any plan- 
etary candidates with host stars that did not obey Equation 
(Q~|) for consistency with the simulated stellar population. This 
step filtered out about half of the sample; 1,135 candidates 
remained. Lastly, we made additional cuts based on planet 
radius, period, and SNR: 

P < 200 days, 
1.5 ^ e <R < 30fl e , (2) 
SNR(^Q8) > 11.5. 

These specific cuts were performed in order to choose a sam- 
ple of planetary candidates with period, radius, and SNR that 
are unlikely to be missed by the Kepler detection pipeline. It 
is necessary to choose a sample with a high degree of com- 
pleteness because we want to compare this sample to simula- 
tions, which are 100% complete. For orbital periods, a range 
extending to 200 days ensures that multiple transits are likely 
to have been observed during Ql-6, which covers approxi- 
mately 486.5 days. For planet radii, the lower bound on R of 
1.5 /?® was chosen because the sample of observed planets 
with r adii R< 1.5 Ra, has much lower d etection efficiencies 
(e.g., IHoward et anl2012t lYoudinll20Tlh . For SNR, we re- 
quired a SNR> 10 for Ql-6, which corresponds to SNR> 1 1 .5 
for Ql-8 by assuming that SNR roughly scales as \^N, where 
N is the number of observed transits. This scaling is per- 
formed because the SNRs of observed KOIs ha ve been re- 
ported for Ql-8, not Ql-6, in lBatalha etaLl ( 120121) . 

Now we discuss how we calculated detection efficiencies, 
and how they were used to debias this filtered sample of plan - 
etary candidates. Following the formalism of lYoudinl d201 ll) . 
for any given planet the net detection efficiency r\ is the prod- 
uct of (1) the geometrical selection effect 7y ge om and (2) the 
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Kepler photometric selection effect ?7phot- The geometrical se- 
lection effect is due to the probability of transit as planets with 
longer periods or larger semi-major axes will have a lower 
probability of crossing in front of the disk of the star as seen 
from Kepler, and is calculated as 77 ge0 m = R*/a- The Kepler 
photometric selection effect is due to the photometric qual- 
ity of the star for detecting planets with specified radius and 
period. Consequently, 77 p h t is a function of both planetary ra- 
dius and period. For a planet with radius R and period P, we 
calculated r^ot as the fraction of stars that can detect such a 
planet with SNR>10, based on the KIC sample obtained in 
the previous section. Lastly, we calculated the net detection 
efficiency as rj = % e om»? P hot- 

We applied the above methods for obtaining detection effi- 
ciencies to our KOI subset of observed planetary candidates. 
For each candidate planet (or detection) in our KOI subset, 
we calculated its net efficiency rj given its observed period P 
and radius R. Since rj is the ratio of the number of detectable 
events to the number of actual planets, the inverse of 77 is equal 
to an estimate of the actual number of planets represented by 
each detection. After we calculated the value of 1 jr\ for each 
detection in our KOI subset, we created debiased, binned his- 
tograms for P and R representing an estimate of the actual 
number of planets per P or R bin. The normalized versions of 
these histograms are given in FigureQ] We used these discrete 
distributions representing debiased planet periods and radii in 
order to randomly select values for each simulated planet's 
period and radius. 

Our methods of obtaining photometric detection effi- 
ciencies (77phn t ) are the same as the methods described in 
iHoward et al.l (120121) used to obtain binned (in P an d R) ef - 
ficiencies that were later analytically fit by lYoudinl d201 ll) . 
except for the following differences. First, we applied these 
methods to Ql-6 data (previous work used Q2 only). Sec- 
ond, we avoided any loss of information that can occur dur- 
ing the binning process and/or the power law fitting process 
by skipping those steps and directly calculating the detection 
efficiency for each individual KOI candidate in our filtered 
sample. Third, we calculated SNR differently by actually us- 
ing each KIC star's observing history during Ql-6 (e.g., not 
all stars are observed each quarter) as well as some other mi- 
nor differences; our m ethods for calculating SNR are detailed 
below in Section IXTl 

2.3. Planet Multiplicity Distribution 

We chose to represent the multiplicity distribution of sim- 
ulated planetary systems-the distribution of the number of 
planets per star-by either a delta-like distribution, a modified 
Poisson distribution, or a bounded uniform distribution. All 
distributions are described by a single parameter. We did not 
consider other distributions such as an expon ential distribu- 
tion, w hich may result in a bad fit to the data dLissauer et al.l 
1201 lbl) . All stars in our model populations are assumed to 
have at least 1 planet, and so our study cannot directly address 
the question of planet occurrence. 

A delta-like distribution is a discrete distribution described 
by a single parameter D. If the parameter D represents an in- 
teger value, than all planetary systems following this distribu- 
tion will have D planets each. If the parameter D does not rep- 
resent an integer value, then each planetary system can have 
either [DJ (floor function) or |~D] (ceiling function) planets. 
For instance, a model population with D = 3 means that all of 
its planetary systems have 3 planets per star, and a model pop- 
ulation with D = 4.75 means that 25% of its planetary systems 
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Figure 1. Observed and debiased distributions of orbital peri- 
ods and radii for a subset (Equation (fJJ) of the planets de- 
tected by Kepler. The top plot for orbital period shows a 
distribution for 0-200 days divided into 60 equal-sized bins. 
The observed distribution shows increasing planet frequency 
at shorter periods, mainly due to the greater geometric prob- 
ability of observing transiting planets at shorter periods. The 
bottom plot for planet radius shows a distribution for 1.5—30 
R® divided into 30 equal-sized bins. 

have 4 planets and 75% of its planetary systems have 5 plan- 
ets. As a result, D also represents the mean number of planets 
per syst em. This is the same as the "uniform" distribution 
used bv lLissauer et all d2011bl) . 

A modified Poisson distribution described by parameter A 
is a discrete distribution that is equivalent to a regular Poisson 
distribution with mean value A, except that we ignore zero val- 
ues. Therefore, by using a modified Poisson distribution we 
restrict the regular Poisson distribution so that values picked 
from the distribution must be greater than (so that each plan- 
etary system must have at least 1 planet). As a result, the mod- 
ified Poisson distribution is not strictly Poisson since its mean 
value is different from A, hence the name "modified Poisson." 
Mathematically, the regular Poisson distribution can be writ- 
ten as 



AV 



k\ 



(3) 



which gives the probability of obtaining a discrete value k 
from such a distribution. 

A bounded uniform distribution is represented by a single 
parameter A. For each planetary system, first a modified Pois- 
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son distribution (as defined above) with parameter A is used 
to draw a value N max (maximum number of planets). Second, 
a discrete uniform distribution with range 1-A^ max is used to 
draw a value representing the number of planets in that sys- 
tem. 

We varied the values of the parameter D for the delta-like 
distribution and the values of the parameter A for the modi- 
fied Poisson and bounded uniform distributions. In total, we 
explored D = 2-l and A = 1-6, each with intervals of 0.25, 
resulting in a total of 63 possibilities for the multiplicity dis- 
tribution. 

2.4. Planet Inclination Distribution 

We chose the inclination distribution in each planetary sys- 
tem to be one of three possibilities: an aligned inclination dis- 
tribution, a Rayleigh distribut ion, or a Rayleigh o f Rayleigh 
distribution, as also used by Liss auer et al.l (1201 lbl) . 

An aligned distribution is the straightforward case where 
relative inclinations in the system are 0°. This is the same as 
a perfectly coplanar system. 

A Rayleigh distribution is a continuous distribution de- 
scribed by a single parameter a, which determines the mean 
and variance of the distribution. Its mathematical form is 



P(k) = 



-k 2 /2a 2 



(4) 



which gives the probability of obtaining a value k from this 
distribution. 

A Rayleigh of Rayleigh distribution, with a single parame- 
ter a a , means that we draw from two Rayleigh distributions 
for each planetary system. First, the given parameter a a de- 
fines the first Rayleigh distribution from which we draw a 
value for a for each planetary system. Second, the drawn 
value of a is then used to define the second Rayleigh distribu- 
tion from which we draw values for the inclinations of planets 
in that system. This allows for the possibility that planetary 
systems in a particular model population have Rayleigh incli- 
nation distributions with different mean and variance. 

We varied a for the Rayleigh distribution and a a for the 
Rayleigh of Rayleigh distribution with the same range of val- 
ues: 1-10° with intervals of 1°, 15°, 20°, and 30°. In total, 
all possibilities for the 3 inclination distributions added up to 
a total of 27 possible inclination distributions. 

2.5. Other Planetary Parameters 

So far we have defined the distributions used in drawing the 
orbital periods, radii, multiplicity, and orbital inclinations of 
simulated planetary systems. Now we discuss all other rel- 
evant planetary parameters: planet mass, orbital semi-major 
axis, eccentricity, and longitude of the ascending node. To 
determine planet mass M, we s caled from its radius R using 
M = R 206 dLissauer et alJl201 lbl) . where M is in units of M ffi 
and R is in units of R®, This commonly-used power law was 
derived by fitting to Earth and Saturn. Each planet's semi- 
major axis is determined using the orbital period and Kepler's 
Third Law. For eccentricity, we assumed that all planets have 
circular orbits. Kepl er candidates do not show evidence for 
large eccentricities (Moorhea d et alJl201 ll) . but we note that 
even moderate eccentricities (<0.25) could affect transit du- 
rations in a manner similar to that caused by differences in 
inclinations of up to a few degrees. For the longitude of the 
ascending node of each planet's orbit, we drew an angle from 
a uniform distribution with values between 0° and 360°. 



Now we have described all relevant stellar and planetary 
parameters to build a fictitious planetary system. In the next 
subsection, we describe stability constraints placed on these 
systems. 

2.6. Stability Requirements 

The first stability constraint is Hill stability, and we describe 
our implementation using the following example. In the con- 
struction of each planetary system, we add a first planet with 
its parameters chosen according to the methods already de- 
scribed. If the system's assigned multiplicity is greater than 
one, we add a second planet. The second planet is accepted 
only if it is Hill stable with the first planet, which occurs when 



fl2 — fli 



R 



> 3.46, 



HL2 



(5) 



where a is the semi-major axis and Rm,2 is the mutual Hill 
radius defined as 



R 



HI, 2 '■ 



Mi+M 2 
3M» 



1/3 



fll +fl2 



(6) 



(e.g.. iGladmanll 1 9931: IChambers et al.lll996l) . Here, M repre- 
sents mass and subscripts 1 and 2 denote the inner and outer 
planets, respectively. If the second planet is not accepted be- 
cause it violates Hill stability with the first planet, then the 
second planet has its radius and orbital period (or equiva- 
lently, mass and semi-major axis) re-drawn until Hill stabil- 
ity is satisfied. We require that all adjacent planet pairs sat- 
isfy this Hill stability requirement. The assumption that ad- 
jacent planet pairs should be Hill stable is reasonable; out of 
the 885 planet candidates located in multi-planet systems de- 
tected during Ql-6 of the Kepler mission, o nly two pairs of 
planet s are found to be apparently unstable by Fa brvckv et alj 
( 120121) on the basis of numerical integrations using nominal 
mass-radius scalings and circ ular orbit assumptions. 

Following the methods of Lissau eFet alj d201 lbl) . we im- 
posed a second stability constraint that any adjacent planet 
trio (any three neighboring planets) must satisfy 



A 



■inner pair of trio 



^outer pair of trio 1 8. 



(7) 



As described for the example above for Hill stability cases, if 
the above constraint in Equation is not satisfied, a planet's 
properties are re-drawn until this stability requirement is ful- 
filled. 

For the creation of multi-planet systems, we added planets 
to each system one at a time, where successful addition of 
a planet occurred when it fulfilled both stability criteria de- 
scribed above. For systems with a large multiplicity, at times 
it was not possible to satisfy stability constraints given our pe- 
riod range and therefore we allowed a maximum of 1000 at- 
tempts per star before discarding it. Therefore, in rare cases, 
model populations have a slightly lower number of planetary 
systems. 

2.7. Summary 

To summarize this section: we have described the steps 
taken to assemble various model populations, each consisting 
of about 10 6 planetary systems. Each of these model popula- 
tions follows different combinations of ( 1 ) multiplicity and (2) 
inclination distributions. Taking all combinations of the val- 
ues for the tunable parameters in these distributions, there are 
a total of 1,701 possibilities. Our goal is to take these model 
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populations, run them through synthetic observations, and see 
which planets are transiting and detectable. The properties of 
these detectable planets will then be compared to the actual 
observed Kepler planets to determine which model produces 
acceptable fits. This allows us to determine the nature of un- 
derlying multiplicities and inclinations. In the next section, 
we discuss how we evaluate these model populations. 

3. EVALUATING MODEL POPULATIONS 

We describe our methods for evaluating each model popu- 
lation: (1) determining detectable, transiting planets, and (2) 
comparing them to the observed Kepler planets using statisti- 
cal tests. 

3.1. Determining Detectable, Transiting Planets 

For each model population, we determined which of the 
planets in each planetary system were transiting. From the 
group of transiting planets, we determined which transiting 
planets were detectable based on SNR requirements. The fol- 
lowing paragraphs describes these methods. 

For each planetary system, we evaluated whether any of the 
planets was transiting. First we picked a random line-of-sight 
(i.e., picking a random point on the celestial sphere) from the 
observer to the system, and the plane normal to this line-of- 
sight vector determines the plane of sky. We then computed 
the planet-s tar distance projected on the plane of sky (e.g., see 
IWinnll20Toh . The minimum in that quantity is the impact pa- 
rameter, which we normalized by the stellar radius. A planet 
transits when its impact parameter is less than 1 . 

We calculated the SNR of transit events in order to deter- 
mine whether or not each transiting planet would have been 
detectable during Ql-6 of the Kepler mission. To do so, first 
we randomly picked a MJD (Modified Julian Date) between 
the first and last cadence mid-times of Ql-6 by using the dates 
provided in the Kepler Data Release Note£| (DRN) 14 and 16. 
We assumed that the planet transited on this randomly-picked 
date, and we explored backwards and forwards in time us- 
ing the planet's orbital period in order to determine the dates 
of all transits during Ql-6. Next, we eliminated any transits 
that occurred in the gaps between quarters; such gaps are nec- 
essary to roll the spacecraft 90° to keep its solar arrays il- 
luminated. In addition, not all stars are observed each quar- 
ter, for instance, due to CCD module failure. As a result, we 
used the planet's host star o bserv ing log that recorded all ob- 
served quarters (see Section |2~TT >, and eliminated any transits 
that occurred in quarters where the host star was not observed. 
Next, we eliminated a small fraction of remaining transits in 
a random fashion to account for a 95% duty cycle because of 
observing downtimes (e.g., breaks for data downlink). After 
accounting for these issues, all remaining N transits for this 
planet were then counted and used to calculate the detection 
SNR as 



(8) 



SNR = | A ^L. 

K*J OCDPP 



where the first fraction represents the depth of the transit and 
ccdpp represents the noise (quadratically interpolated to the 
transit duration length using the set of measured {3-hr, 6-hr, 
12-h r} CDPP values) for the planet's host star (see Section 
|2~T1 i. If the SNR is at least the value of the SNR threshold for 

2 The DRN are available for download from MAST at 
|http: / /archive . stsci . edu/kepler/data_release.html| 



detection (SNR>10), then the transiting planet is defined as 
detectable. 

3.2. Comparing to Observations 

Here we describe how we evaluated each model popula- 
tion's goodness-of-fit when compared to observations. All 
model populations are distinct from each other due to the 
different parameter values chosen for their underlying mul- 
tiplicity and inclination distributions. To evaluate how well 
each model population matches the observations, we com- 
pared them using 2 observables: (1) observed multiplicity 
vector and (2) observed distribution of normalized, transit du- 
ration ratios. These observables were computed for both the 
simulated planets and the KOI candidates. 

We discuss these 2 types of observables in greater detail. 
The multiplicity vector fi represents the number of systems 
observed with j transiting planets, where j = 1,2,3,4,5,6,7+. 
In other words, the multiplicity vector describes how many 
systems are observed to have a single planet, how many sys- 
tems are observed with 2 planets, etc. Systems with 7 or more 
planets are placed in the same multiplicity category: j = 7+. 
The normalized, tran sit duration ratio £ is defined as (e.g., 
Fabrvcky 'eTani2012h 



Tdur.l/f] 



1/3 



1/3 : 



(9) 



where the transit duration Ta m is normalized by the orbital pe- 
riod raised to the 1/3 power P 1 ^. Subscripts 1 and 2 represent 
any pair of planets in the same system where 1 denotes the 
inner planet (the planet in the pair with the shorter orbital pe- 
riod) and 2 denotes the outer planet. A /'-planet system has 
a total of l)/2 pairs. We decided to fit to values of £ 
because they contain information on mutual inclinations of 
observed multi-planet systems. For instance, coplanar sys- 
tems tend to have values of £ > 1, because the inner planet's 
impact parameter i s smaller and hence it s normalized transit 
duration is longer. Fabr vckv et al.l d2012l) used the observed 
Kepler distribution of £ to determine that the observed mutual 
inclinations are low at 1.0-2.3°. 

We calculated the (p,, Q obs erva bles for the subset of KOIs 
previously obtained in Section l2T2l We also calculated the (pi, 
£) observables for each model population by considering all 
detectable, transiting planets. 

Our next step was to obtain a statistical measure of each 
model population's goodness-of-fit to the observed popula- 
tion by calculating a test statistic and its corresponding prob- 
ability. A larger probability indicates a better match to ob- 
servations. To compare multiplicity vectors fj,, we performed 
chi-square tests, which are appropriate for discrete data. Chi- 
square tests are valid when the chi-square probability distri- 
bution is a good approximation to the distribution of the chi- 
square statistic (Equation ( fT0T > below). As shown by theoret- 
ical investigations, the approximation is usually good when 
there are at least 5 discrete categories and at least 5 values 
in ea ch category for the model population's /x (e.g., iHoell 
11984ft . These requirements may not always be satisfied in 
our case, and so we also perform a randomization method 
to calculate chi-square probabilities, as described below. To 
compare £ distributions, we performed Kolmogorov-Smirnov 
(K-S) tests, which can be used for continuous data. K-S tests 
have no restrictions on sample sizes, and are also distribution- 
free and so therefore can be applied for any kind of underly- 
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ing di stribution of the comparison samples (e.g. JStuart & Ordl 
1991). These tests are described in more detail below. 

The chi-square test used for comparing multiplicity vectors 
uses a chi-square statistic x 2 computed as (e.g., iPress et all 
fl992h 

J 1 

where Oj represents an observed quantity, Ej represents an 
expected or theoretical value, and j represents the index of the 
multiplicity vector (i.e., 1 for 1-planet systems, N for N-planet 
systems, etc) being summed over. Scaling was performed to 
adjust the model multiplicity vector such that Oj = ^Ej. 
Large values of x 2 represent worse fits and larger deviations 
between Oj and Ej quantities, meaning it is not likely the Oj 
values are drawn from the population represented by the Ej 
values. For comparison between each model population's ex- 
pected values and observed values, we only considered j in- 
dices where Ej ^0 (noting the appearance of Ej in the denom- 
inator of Equation (TTOb>. Since we assigned no weight to cat- 
egories with Ej = 0, we note that there may be cases (i.e., for 
the delta-like multiplicity distribution) where a given model 
population with only 1-planet and 2-planet systems may yield 
a low x 2 even though it is in fact not a good match to the 
data because it does not predict any 3+ planet systems that 
are actually observed. From our results, we find that model 
populations based on delta-like multiplicity distributions do 
not provide a good match to the data and are therefore not 
discussed further; consequently none of the x 2 probabilities 
presented in this paper suffer from this issue. 

We computed two chi-square probabilities for each calcu- 
lation of the statistic in Equation (TTOb . First, the standard chi- 
square probability is an incomplete gamma function that can 
be computed with knowledge of x 2 and the number of degrees 
of freedom. For our case, the number of degrees of freedom 
is equal to the number of indices with Ej / subtracted by 1 
(to account for the scaling constraint). This standard proba- 
bility function is a good approximation as long as the number 
of bins or the number of values in each bin is large. Given that 
some higher- j categories in our multiplicity vector may have 
low numbers, we performed a second calculation of the chi- 
square probability using a more robust randomization method 
(e.g., [Good 2006). In this method, we determined how often 
we could pick a new distribution from the model distribution 
with xi evj mo dei (comparing new and model distributions) that 
was worse than the x^ bs model (comparing observed and model 
distributions). In practice, we randomly picked 10,000 times 
per model population and the fraction of them that yielded 
Xnew,modei > Xobs.modei represented the chi-square probability. 
We found that in most cases, there was reasonable agree- 
ment between the chi-square probabilities from the standard 
function and from the randomization method. The chi-square 
probabilities quoted later in this paper are those obtained from 
the randomization method. 

For the K-S test used for comparing £ distributions, we 
calcul ated the K-S st atistic and the standard K-S probability 
(e.g.. lPress et al.lfl992l) . The K-S statistic D is the maximum 
difference of the absolute value of two distributions given as 
cumulative distribution functions, or 

D = max\0 Nl (0-E N2 (0\, (11) 



where Ojv, (0 represents the observed cumulative distribution 
function of £ values, represents the expected or model 

cumulative distribution function, and N\ and N2 is the num- 
ber of £ values for the observed and expected distributions, 
respectively. The significance or probability of D represents 
the chance of obtaining a higher value of D, and can be cal- 
culated with knowledge of D and the number of points N\ 
and N2- This test is valid even with an unequal number of 
points, i.e., N\ / N2. We note that N2 can vary significantly 
between model populations and this can affect the reliability 
of the calculated probability, such as when comparing a prob- 
ability obtained from one model population to the probability 
obtained from another model. This is most pertinent for cases 
with large differences in N2 (i.e., between model populations 
with small and large multiplicities), where a poor fit is harder 
to discern when there are fewer model values N2 to compare 
to observed values. This issue can manifest itself as outliers in 
our results for model populations with fewer number of plan- 
ets per system. As we will show, we do not find any evidence 
of such outliers in our results. 

In the next section, we describe our calculations of these 
probabilities to evaluate the match between each model pop- 
ulation and the data. 

4. RESULTS 

In this section, we present the results for the fit between 
each model population and the KOI subset that obeys Equa- 
tions ([TJ and ©, and these results are illustrated in Figures 
|2]-|4] We do not show any results incorporating the delta-like 
multiplicity distribution because its results do not provide a 
good match to the data. We will describe each of these figures 
in turn. 

Figure [2] shows the probabilities resulting from the x 2 test 
that compared the observed and model multiplicity vectors 
(number of 1-planet systems, number of 2-planet systems, 
etc.). Most models can be ruled out with 3-sigma confidence, 
and these tend to be systems with many planets and low in- 
clinations or systems with few planets and high inclinations. 
On the other hand, models that are consistent with the data 
include systems of few planets with low inclinations or sys- 
tems of many planets with high inclinations; this degeneracy 
is evident in the upward slopes in Figure [2] Models with a 
bounded uniform multiplicity distribution and a Rayleigh of 
Rayleigh inclination distribution appear to provide a better fit 
(higher x 2 probabilities) than other models. 

Figure [3] depicts the probabilities from the K-S test that 
evaluated the fit between the observed and model £ (Equa- 
tion (0) distributions of multi-planet systems. The distribu- 
tion of £ values is used to gauge the extent of coplanarity or 
non-coplanarity in multi-planet systems. This diagram shows 
that higher probabilities are obtained for models with low- 
inclination distributions. Nearly all models with a Rayleigh 
a or Rayleigh of Rayleigh a a parameter greater than 3° can 
be ruled out with 3-sigma confidence. None of the models 
with perfectly aligned (coplanar) systems provide acceptable 
fits. The best fits with highest probabilities are Rayleigh or 
Rayleigh of Rayleigh inclination distributions with a, a a = 1 °. 
Comparison between Figures|2]and[3]shows that the degener- 
acy evident in Figure [2] is broken when considering Figure [3] 
as well-by also examining Figure[3] we see that only systems 
with fewer planets and lower inclinations are consistent with 
the data. 

Figure|4]illustrates the combined probabilities from Figures 
[2] and [3] The combined probability is the product of the x 2 



6 



X 2 Probabilities for Multiplicity 




2.25 

2 

1.75 
1.5 
1.25 

1 | 

1 234567891015 20 30 1 234567891015 20 30 
Al. I Ra. I I RR. I 

Underlying Inclination Distribution 

Figure 2. Diagram showing probabilities from x 2 tests that evaluated the goodness-of-fit of various models with respect to planet 
multiplicity. Each square on this grid represents a different model, and each square's color depicts the probability (the higher the 
probability, the better the fit). The x-axis shows the underlying inclination distribution, and can be aligned (AL), Rayleigh (Ra.), 
or Rayleigh of Rayleigh (RR.). For Al., the number given is the inclination, or 0°. For Ra. and RR., the numbers indicate the a 
or a a parameter (in degrees) for the distribution. The y-axis shows the underlying multiplicity distribution, and can be modified 
poisson (Po.) or bounded uniform (BU.). For Po, the number given is A of the Poisson distribution and represents approximately 
the mean number of planets per system. For BU., the number given is A of the Poisson distribution from which the maximum 
number of planets is drawn. This plot can be approximately divided into 4 quadrants, where in each quadrant the y-axis has the 
number of planets per system increasing to the top and the x-axis has inclinations increasing to the right. 
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K-S Probabilities for £ 
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Underlying Inclination Distribution 

Figure 3. Same as Figure |2j except this diagram shows K-S probabilities of £ distributions. These probabilities represent the 
goodness-of-fit of various models with respect to orbital inclinations. 
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Combined Probabilities 
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Underlying Inclination Distribution 

Figure 4. Same as Figure [2] except this diagram shows combined (multiplied) probabilities from x 2 an d K-S tests. These proba- 
bilities represent the overall goodness-of-fit of various models to the Kepler sample. 
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and K-S probabilities. We make the assumption that these 
two probabilities are independent. This product of two prob- 
abilities is shown in Figure [4] from which it is evident that 
most models can be ruled out as unacceptable fits. Degenera- 
cies present in both Figures [2] and [3] are broken when the data 
in these figures are combined in Figure [4] The best fits are 
models with lower inclinations and lower multiplicities, and 
in particular, the fits with a bounded uniform distribution in 
multiplicity and a Rayleigh or Rayleigh of Rayleigh distribu- 
tion in inclination are most consistent with the data. 

4.1. Overall Best-Fit Models: 
Few Planets and Low Inclinations 

As seen in Figure [4] the overall best-fit models with 
the largest combined probabilities are models with low- 
multiplicity bound uniform distributions (represented by A) 
and with low-inclination Rayleigh (represented by a) or 
Rayleigh of Rayleigh (represented by a a ) distributions. The 
best-fit model with a bounded uniform and Rayleigh distri- 
bution has A = 2.50 and a = 1°. The best- fit model with a 
bounded uniform and Rayleigh of Rayleigh distribution has 
A = 2.25 and a a = 1°. There are no qualitative differences be- 
tween the quality of these two fits. Based on these two best-fit 
models, we find that 75-80% of planetary systems have 1 or 

2 planets with orbital periods less than 200 days. In addi- 
tion, over 85% of planets have orbital inclinations less than 

3 degrees (relative to a common reference plane). These re- 
sults represent our best estimate of the underlying distribu- 
tions of Kepler transiting systems for the parameter space ex- 
plored here (Equations (UJ and (0). Assuming that the Kepler 
sample is representative, these distributions also describe the 
architecture of planetary systems in general. 

We analyze the best-fit model with a bounded uniform 
(A = 2.25) and a Rayleigh of Rayleigh (a a = 1°) distribution 
in greater detail. First, in Figure [5] we plot the underlying 
multiplicity and inclination distributions. The top and middle 
panels of the figure show the first and second steps of picking 
multiplicity and inclination values from the distributions. The 
bottom plots shows the resulting histograms of values picked 
from these distributions. Second, we show the comparisons 
between the best-fit model and the data. Table[T]compares the 
numbers of observed transiting systems to the numbers of de- 
tected transits from the best-fit model. There is a reasonable 
match in each multiplicity index and an overall match proba- 
bility of 25.2%. Figure [6] illustrates the comparison between 
the distribution of £ values from Kepler observations and the 
distribution of £ values from simulated detections using the 
best-fit model. The numbers in each bin have been normal- 
ized to enable comparison, and the match between distribu- 
tions has a probability of 27.5%. Taken together, the prob- 
abilities from each fitted observable (\ 2 probability for the 
multiplicity fit and K-S probability for the £ distribution fit) 
are multiplied to yield a combined probability of 6.9%. 

The low inclinations in our best-fit models have interest- 
ing implications for planet formation and evolution theories. 
Our findings that these systems are relatively coplanar (at least 
out to 200 days) are in favor of standard models that suggest 
planet formation within a protoplanetary disk. In addition, 
strong influences by external perturbers such as Kozai pro- 
cesses, scattering, or resonances are not likely to play a ma- 
jor, lasting role given that these systems do not generally have 
high, excited inclinations. Theories of planet formation and 
evolution can be tested against the architecture of planetary 



Table 1 

Multiplicity Vector: 
Number of Systems with j Detectable Planets 



Name j = 


1 


2 


3 


4 


5 


6 


7+ 


Observed 


542 


85 


24 


4 


1 


1 





Best-fit 


540.6 


92.9 


19.1 


3.6 


0.6 


0.2 


0.0 



Comparison of the multiplicity vector between observations and 
results from one of t wo b est-fit models, here with A = 2.25 and 
(T a = 1° (see Section |4~TY The \ 2 probability for this match is 



systems "illuminated by the Kepler mission. 

4.2. Effects of False Positives and Survey Incompleteness 

The results of our study do not take into account the ex- 
istence of false positives that mimic real planet transit sig- 
nal s in the KOI dataset. For the Kepler sample re leased 
by iBorucki et ail (1201 lbh . iMorton & Johnson! (1201 ll) esti- 
mated that false positive probabilities were gen erally low at < 
10%, a nd an analysis using binary statistics by Lissaue r et al.l 
(120121) estimated that ^98% of the multi-plan et systems are 
real. For the latest Kepler sa mple released by Batal ha et al.l 
(120121) . iFabrvckv et al.l (120121) discussed how ~96% of pairs 
in multi-transiting candidate systems are likely to be real plan- 
ets using stability arguments. For the purposes of our study, 
we assumed that all planetary candidates are real. We expect 
that the existence of false positives can affect the fit between 
multiplicity vectors (e.g., Tabled, particularly for the case of 
singly transiting systems where the actual population may ac- 
tually have lower numbers than reported due to false positives. 
We do not expect that false positives will significantly influ- 
ence the fit between normalized transit duration ratios (e.g., 
Figure[6]i since such ratios are only computed for multi-planet 
systems, which we expect to be a higher-fidelity sample than 
singly transiting systems. 

Another effect that can impact our results is Kepler sur- 
vey incompleteness. The Kepler team is currently investi- 
gating the completeness of the pipeline by inserting artifi- 
cial tr ansits into the pipelin e and determining their recovery 
rates. Batalha et al. (2012) suggested that the KOI catalog 
released in 2011 by Borucki et al. (2011g) may suffer from 
incompleteness issues because the planet gains seen between 
the 2011 and 2012 catalogs cannot be completely explained 
by the longer observation window, and that it is possible the 
2012 catalog may still suffer from some incompleteness is- 
sues. This can affect our results in terms of the comparison 
between the observations and our model's computed observ- 
ables, and we expect that its impact is a function of SNR and 
other parameters. Its impact will be better known after the 
Kepler team's detailed study of pipeline completeness has fin- 
ished, and pipeline completeness is expected to improve in the 
future with planned pipe line upgrades that are already under- 
way (Batalha et al. 2012). 

5. COMPARISON WITH THE SOLAR SYSTEM 
AND WITH PREVIOUS WORK 

We compare the results of our study with the properties of 
planets in the Solar System as well as with some relevant pre- 
vious studies. 

5.1. Results Are Consistent With the Solar System 

Based on the underlying multiplicity distribution of our 
best-fit models, 75-80% of planetary systems have 1 or 2 plan- 
ets with orbital periods less than 200 days (see Section |4~TT i. 
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Underlying Multiplicity Distribution 

Step 1 




Underlying Inclination Distribution 

Step 1 




Step 2 



Step 2 




Number of Planets per System 



Result 




4 6 8 

Number of Planets per System 



Result (Cumulative Plot) 



^0.9 




Number of Planets per System 




Planet Inclination [deg] 



Result 



0.08 



>.0.06 



pO.04 



0.02 




Planet Inclination [deg] 
Result (Cumulative Plot) 




Planet Inclination [deg] 



Figure 5. One of two best-fit models for multiplicity {left column) and inclination {right column). The multiplicity distribution is 
a bounded uniform distribution: for each planetary system (1) we draw a value N max from a modified Poisson distribution with 
A = 2.25, and (2) we choose the number of planets by uniformly picking a value between 1 and N max . The resulting distributions 
of multiplicities are shown in the bottom plots. The inclination distribution is a Rayleigh of Rayleigh distribution: (1) for each 
planetary system we draw a value for a from a Rayleigh distribution with parameter a a = 1°, and (2) for each of its planets we 
draw a value for inclination from a Rayleigh distribution with parameter a. The resulting distributions of inclinations are shown 
in the bottom plots. 
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Figure 6. Comparison of the £ distribution between observa- 
tions and results from one of tw o bes t-fit models, here with 
A = 2.25 and <j a = 1° (see Section l4TTT i. The quantity £ is sen- 
sitive to orbital inclinations in multi-planet transiting systems. 
The K-S probability for this match is 27.5%. 

Our analysis is agnostic about the number of planets with 
larger orbital periods or different radius/mass regimes than 
considered here; recall that these limits are due to the param- 
eter space constraints imposed by Equations ((TJ and |0. In 
comparison, the Solar System has 1 planet, Mercury, with an 
orbital period less than 200 days. As a result, our best-fit mul- 
tiplicity distributions are consistent with the Solar System, but 
more data and extension of this study to longer orbital periods 
are warranted before any definitive multiplicity comparison 
can be made with the Solar System. 

Based on the underlying inclination distributions of our 
best-fit models, over 85% of pl anet s have orbital inclinations 
less than 3 degrees (see Section [4Tt , suggesting a high degree 
of coplanarity. This is compatible with the inclinations seen 
among the planets in the Solar System, if we allow ourselves 
to extrapolate beyond the period limits of our study. In the 
Solar System, 7 out of the eight planets (or 87.5%) have in- 
clinations less than 3 degrees relative to the invariable plane, 
with Mercury as the exception. 

5.2. Comparison with Previous Studies 

iLissauer et al.l (1201 lbl) investigated the dynamical proper- 
ties of_jnvd^i-2]^ets^stems in the KOI catalog announced 
by iBorucki et alJ (1201 lbl) . They also presented a detailed 
analysis of the inclinations of Kepler planetary systems. To 
accomplish this, they created a host of planetary models 
following different planet multiplicity and inclination distri- 
butions, determined which of those planets were transiting 
and detectable, and compared the resulting transit numbers 
with the observed transit numbers to determine each model's 
goodness-of-fit. Differences between their methods and ours 
include different methods for obtaining debiased period and 
radius distributions, different radius and period parameter 
ranges, an observed £ distribution that was used as a fitted 
observable in our study, and different d ata sets (we used th e 
most recently rele ased KOI data set: | Batalha et all 120121) . 
From their results, Lissa uer et al.l (1201 lbl) ruled out systems 
with small numbers of planets and large inclinations as well as 
systems with large numbers of planets and small inclinations, 
which we also found to be inconsistent with the data (see Fig- 
ure |2j. In addition, they found a degeneracy in their results 
between underlying distributions of inclination and multiplic- 



ity; the number of planets per system could be increased if the 
inclination was also increased and still provide good fits to the 
data, also seen in our Figure|2] Ultimately they discounted the 
case of thick systems with many planets by invoking results 
from radial velocity surveys that s uggested low inclinations 
(see discussion in Fabr ycky et al.ll2012h . In our study, we 
were able to reject that scenario because it provided poor fits 
to the observed £ dist ri bution (see Figure [3]). 

Tre maine & Dongl (120121) applied a general statistical 
model to Kepler as well as to radial velocity surveys to an- 
alyze the multiplicity and inclination distributions. They con- 
cluded that Kepler data alone could not place constraints on 
inclinations, and acceptable possibilities ranged from thin to 
even spherical systems. Systems with large rms inclinations 
could be fit by the data as long as some of them had a large 
number of planets. However, when jointly analyzing both Ke- 
pler and radial velocity data, they were able to place con- 
straints on inclinations to show that mean planetary inclina- 
tions are in the range 0-5°. These relatively low inclina- 
tions are consistent with our results showing how models with 
lo wer inclination s provi ded better fits to the data. 

Figueira et al.1 (120121) determined underlying inclination 
distributions by applying information from both HARPS 
(High Accuracy Radial Velocity Planet Searcher) and Kepler 
surveys. They assumed that if the HARPS and Kepler sur- 
veys share the same underlying population, then the different 
detection sensitivities of these two surveys with their detected 
samples of planets should allow them to determine underlying 
inclinations. To do so, they created a synthetic population of 
planets with a given multiplicity distribution as previously de- 
termined using HARPS data, and planets were given various 
inclination distributions (aligned, Rayleigh, and isotropic). 
Each model's goodness-of-fit was determined by comparing 
the frequency of transiting systems with the Kepler sample. 
They found that the best fits were obtained using models with 
inclinations prescribed by a Rayleigh distribution with a < 1 ° . 
In our stu dy, we used a differe nt data set based on the latest 
release by Batalhaet al. (2012j), we fit for an additional pa- 
rameter (the intrinsic multiplicity distribution), and we fit to 
an extra observable (the resulting £ distribution). Given that, 
our overall results are in agreement since we also found that 
low inclinations with a Rayleigh parameter cr ~1° are consis- 
te nt with the data . 

Fabrycky et alJ d2012l) presented important properties of 
multi-planet candid ates in the KOI catalog released by 
iBatalha et alJ (120121) . They found that almost all of these sys- 
tems are apparently stable (using nominal mass-radius and 
circular orbit assumptions), which reinforces the high-fidelity 
nature of these multi-planet detections. In addition, they de- 
rived the mutual inclinations of observed planets using the 
distribution of £ from the Kepler sample, which we used in 
our study as a fitted observable. They found that mutual in- 
clinations are constrained to the range of 1-2.3° for observed 
planets, and concluded that planetary systems are typically 
flat. Our analysis is different and has some advantages in that 
we proposed underlying planet multiplicity and inclination 
models, ran them through synthetic observations, and fitted 
the models' results to observed transit number s and £ distribu- 
tion. O ur study agrees with the conclusions of Fabrycky et al. 
(120121) that planetary systems tend to be relatively coplanar to 
w ithin a few degrees. 

Joha nsen et al.l d2012l) used a different approach than some 
previous authors. They created synthetic triple planet systems 
by assuming they could be based off of observed Kepler triply 
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transiting systems. They determined which of their synthetic 
planets could be transiting and observable, and then calcu- 
lated the resulting number of singly, doubly, and triply tran- 
siting systems to compare with the actual number of transits 
in the Kepler data. They also repeated these steps for different 
inclinations. From fitting to observed doubly and triply tran- 
siting systems, they found low mutual inclinations of < 5°. 
Although our methods are different, such low mutual inclina- 
tions are compatible with our findings that these systems are 
intrinsically thin. 

6. CONCLUSIONS 

We have investigated the underlying distributions of multi- 
plicity and inclination of planetary systems by using the sam- 
ple of planet candidates discovered by the Kepler mission 
during Quarters 1 through 6. Our study included solar-like 
stars and planets with orbital periods less than 200 days and 
with radii of 1.5-30 R©. We created model populations rep- 
resented by a total of two tunable parameters, and we fitted 
these models to observed numbers of transiting systems and 
to normalized transit duration ratios. We did not include any 
constraints from radial velocity surveys. Below we list the 
main conclusions of our study. 

1. From our best-fit models, 75-80% of planetary systems 
have 1 or 2 planets with orbital periods less than 200 days. 
This represents the unbiased, underlying number of planets 
per system. 

2. From our best-fit models, over 85% of planets have or- 
bital inclinations less than 3 degrees (relative to a common 
reference plane), implying a high degree of coplanarity. 

3. Compared to previous work, our results do not suffer 
from degeneracies between multiplicity and inclination. We 
break the degeneracy by jointly considering two types of ob- 
servables that contain information on both number of planets 
and inclinations. 

4. The underlying multiplicity distribution is consistent 
with the number of planets in the Solar System with orbital 
periods less than 200 days. If we extrapolate to beyond 200 
days, we find that the underlying distribution of inclinations 
derived here is compatible with inclinations in the Solar Sys- 
tem. Our results are also consistent with the standard model of 
planet formation in a disk followed by evolution without any 
strong perturbative force that could have caused a significant 
and lasting impact on inclinations. 

Continued observations by the Kepler mission will improve 
the detectability of new candidate planets covering a larger 
swath of parameter space, especially to longer orbital periods 
and smaller planetary radii. We anticipate that future statisti- 
cal work will further boost our understanding of the underly- 
ing architecture of planetary systems. 



We thank Dan Fabrycky for useful discussions, as well as 
the entire Kepler team for procuring such an excellent dataset 
of planetary systems. 
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